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We study the statistical distributions of the spins of generic black-hole binaries during the inspiral 
and merger, as well as the distributions of the remnant mass, spin, and recoil velocity. For the 
inspiral regime, we start with a random uniform distribution of spin directions Si and 5*2 over the 
sphere and mag nitudes \Si/ml\ = |52/mi| = 0.97 for different mass ratios, where Si and rui are the 
spin-angular momentum and mass of the ith black hole. Starting from a fiducial initial separation 
of Ti = 50M, we perform 3.5-post-Newtonian-order evolutions down to a separation of r/ — 5M, 
where M = mi -|- m2, the total mass of the system. At this final fiducial separation, we compute 
the angular distribution of the spins with respect to the final orbital angular momentum, L. We 
perform 16* = 65536 simulations for six mass ratios between q = 1 and q = 1/16 and compute 

the distribution of the angles L ■ A and L ■ S, directly related to recoil velocities and total angular 
momentum. We find a small but statistically significant bias of the distribution towards counter- 
alignment of both scalar products. A post-Newtonian analysis shows that radiation-reaction-driven 
dissipative effects on the orbital angular momentum lead to this bias. To study the merger of 
black-hole binaries, we turn to full numerical techniques. In order to make use of the numerous 
simulations now available in the literature, we introduce empirical formulae to describe the final 
remnant black hole mass, spin, and recoil velocity for merging black-hole binaries with arbitrary 
mass ratios and spins. Our formulae are based on the post-Newtonian scaling, to model the plunge 
phase, with amplitude parameters chosen by a least-squares fit of recently available fully nonlinear 
numerical simulations, supplemented by inspiral losses from infinity to the ISCO. We then evaluate 
those formulae for randomly chosen directions of the individual spins and magnitudes as well as the 
binary's mass ratio. The number of evaluations has been chosen such that there are 10 configurations 
per each dimension of this parameter space, i.e. 10^. We found that the magnitude of the recoil 
velocity distribution decays exponentially as P{v) ~ exp(— u/2500 km s~^) with mean velocity 
< v >= 630 km s~^ and standard deviation V< > — < f >^ = 534 km s~^, leading to a 23% 
probability of recoils larger than 1000 km s~^, and a highly peaked angular distribution along the 
final orbital axis. The studies of the distribution of the final black-hole spin magnitude show a 
universal distribution highly peaked at Sf/m^f = 0.73 and a 25° misalignment with respect to 
the final orbital angular momentum, just prior to full merger of the holes. We also compute the 
statistical dependence of the magnitude of the recoil velocity with respect to the ejection angle. 
The spin and recoil velocity distributions are also displayed as a function of the mass ratio. We 
finally also compute the effects of the observer orientation with respect to the recoil velocity vector 
to take into account the probabilities to measure a given redshifted (or blueshifted) radial velocity 
of accretion disks with respect to host galaxies. 

PACS numbers: 04.25.Dm, 04.25.Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION 

Astrophysical black-hole (BH) binaries are character- 
ized by the mass ratio q = mi/m2 < 1 of the smaller to 
larger BH, where is the mass of BH i, the total mass 
M = mi -1-7712, eccentricity e (assumed to be very small), 
and spins Si (where Si/mf < 1). In addition, it is often 
convenient to parameterize the binary with the symmet- 
ric mass ratio i] = q/{l + q)^ , specific spins di = Si/mf, 
total spin S = Si + S2, A ~ M{S2/m2 — Si/nii), and 
orbital angular momentum L. 

The relevance of the spins to the dynamics of black- 
hole (BH) mergers was recognized soon after the break- 
through in numerical relativity 1 3 that allowed for the 
long term, stable numerical evolution of such systems. 
Notable examples of the early findings are the 'hangup' 



effect a repulsive spin-orbit interaction, that delays 
the merger of black-hole binaries (BHB) when the spins 
are aligned with the orbital angular momentum, and 
simultaneously causes the system to radiate excess an- 
gular momentum, leading to a remnant BH with sub- 
maximal spin. The same mechanism produces an addi- 
tional attractive effect when the spins are counter-aligned 
with the orbital angular momentum, leading to a prompt 
merger. Thus the radiation of angular momentum and 
energy is asymmetric with respect to the relative orien- 
tations of the total spin angular momentum vector and 
the orbital angular momentum. 

When the spins are not exactly aligned or counter 
aligned, new effects appear (the hangup effect is still 
present). Precession of the spins is important dynami- 
cally because it cause the orbital plane to strongly pre- 
cess just prior to merger The final spin of the merged 
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hole can flip with respect to the directions of the individ- 
ual ones, mainly due to the addition of the orbital angu- 
lar momentum [5]. While the spin-orbit coupling leads 
to strong precessional effects near merger, the magni- 
tudes of the spins are not affected to the same degree. 
In particular spin-orbit interactions are too weak to in- 
duce the binary to corotate (or maintain corotation of 
an initially corotating binary) at the last stages of the 
merger because the timescale for the radiation driven 
inspiral is much smaller than the spin-orbit interaction 
timescale [Bj. 

Numerous other papers have studied different spin ef- 
fects, such as the large recoil velocities acquired by the 
remnant of the merger of two spinning black holes [7Ml2| 
and long term evolutions of generic BHBs (i.e. unequal 
mass and unequal, randomly-oriented spins) [131 to 
cite a few of the nearly one hundred papers published on 
the subject since 2006. 

The characterization of the remnant black hole (BH) 
as the by-product of a generic BH binary (BHB) merger 
is of great astrophysical interest as it allows one to model 
the growth of BHs during the evolution of the universe 
and their effect on the dynamical evolution of galactic 
cores and globular clusters, as well as the collisions of 
galaxies and stellar size binary systems. Thanks to the 
recent breakthroughs in Numerical Relativity [IH3] one 
can now precisely compute the masses, spins and recoil 
velocities of these merged BHBs from fully nonlinear nu- 
merical simulations. 

The modeling of the remnant black hole using fully- 
numerical techniques was pioneered by the 'Lazarus 
method' [15, for spinning black holes followed by the 
breakthrough 'moving puncture' approach. In Refs. [31- 
[B] the authors studied BHBs characterized by equal- 
mass, equal-spin individual BHs, with the spins aligned 
or counter-aligned with the orbital angular momentum, 
using fully nonlinear numerical calculations and found 
a simple ad hoc expression relating the final mass and 
spin of the remnant with the spins of the individual BHs. 
This scenario was later revisited in [TBI Hi] and the for- 
mula for the remnant spin was generalized (by assuming 
that the angular momentum is only radiated along the 
orbital axis, and neglecting the energy loss) in ^3] for 
arbitrary BH configurations (although in the latest pa- 
per of this sequel this condition was removed |19|.) In 
[20j a more general ad hoc fitting function was proposed. 
A more comprehensive approach was proposed in |21j : 
where a generic Taylor expansion, reduced by the phys- 
ical symmetries of the problem, was used to fit the ex- 
isting full numerical simulations. A different approach 
was presented in |22j where the particle limit approxima- 
tion was extended to the equal-mass case and the effects 
of post-ISCO (Innermost Stable Circular Orbit) gravi- 
tational radiation were neglected. This approach was 
further improved in |23j by taking binding energies into 
account. All of these approaches show a certain degree of 
agreement with the remnant masses and spins obtained 
in the few dozen fully nonlinear numerical simulations 



available, but there remains significant uncertainties con- 
cerning their accuracy outside this range of parameters. 
In this paper we propose a set of formulae that incorpo- 
rate the benefits of both approaches in a unified way. 

Due to the large astrophysical interest of computing 
remnant recoil velocities, the modeling of recoil velocities 
followed an independent path, particularly since the dis- 
covery [3 [23] that the spins of the black holes play a cru- 
cial role in producing recoils of up to 4000 km s^^. The 
importance of modeling the recoil velocities as a function 
of the astrophysical parameters of the progenitor binary 
was quickly realized [7| |H| [H] • 

The news that the merger of binary black holes can 
produce recoil velocities up to 4000 km s"^, and hence 
allow the remnant to escape from major galaxies, led 
to numerous theoretical and observational efforts to find 
traces of this phenomenon. Several studies made pre- 
dictions of specific observational features of recoiling su- 
permassive black holes in the cores of galaxies in the 
electromagnetic spectrum [^5H5T] from infrared [32] to 
X-rays [33'-'3S] and morphological aspects of the galaxy 
cores 36 .38j. Notably, there began to appear observa- 
tions indicating the possibility of detection of such ef- 
fects [3PJ?T1, and although alternative explanations are 
possible 42 - 45J , there is still the exciting possibility that 
these observations can lead to the first confirmation of a 
prediction of General Relativity in the highly-dynamical, 
strong-field regime. 

In our approach to the recoil problem [Tj [24] we chose 
to use post-Newtonian theory as a guide to model the 
recoil dependence on the physical parameters of the pro- 
genitor BHB (See Eqs. (3.31) in [IB]), while arguing that 
only full numerical simulations can produce the correct 
amplitude of the effect. Bearing this in mind, we pro- 
posed an empirical formula for the total recoil veloci- 
ties (see Eq. (22) below.) Our heuristic formula describ- 
ing the recoil velocity of a black-hole binary remnant 
as a function of the parameters of the individual holes 
has been theoretically verified in several ways. In [23] 
the cos G dependence was established and was confirmed 
in j47j for binaries with different initial separations. In 
[15] the decomposition into spin components perpendic- 
ular and parallel to the orbital plane was verified, and 
in [49 it was found that the quadratic-in-spin corrections 
to the in-plane recoil velocity are less than 20 km s~^. 
Recently in |50| we confirmed the leading 77^ (where 
Tj = (^™:]_™^^)2 is the symmetric mass ratio) dependence 
of the large recoils out of the orbital plane. 

Since the magnitude (and direction) of the recoil ve- 
locity of the remnant black holes depend so sensitively 
on the spin orientation just around the time of the for- 
mation of a common event horizon, it is important to 
establish that random oriented spins of individual black 
holes at large separations (as a plausible initial astro- 
physical scenario) lead to randomly oriented black holes 
near merger, or if there is some bias in their orientations 
by the time they get very close together (i.e. at typical 
numerical simulations separations of a few M, where AI 
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is the total mass). It has been argued recently [5TH55] 
that the presence of gas and accretion of the individ- 
ual black holes during the inspiral phase for long time 
scales can lead to a preferential alignment of spins with 
the orbital angular momentum, and hence to a configura- 
tion that leads to modest recoil velocities (a few hundred 
km s^^). The latest results by Dotti et al [S3] indicate 
that spin alignment occur in the scale of a few million 
years within 10 degrees of the orbital angular momen- 
tum for cold disks and 30 degrees for warmer disks. The 
proportion of wet to dry mergers in the universe still 
needs to be established. Current rough estimates give 
comparable percentages for both kinds of mergers. With 
3 to 5 mergers per galaxy during their lifetime for current 
spiral and elliptic galaxies respectively |102j . While the 
verification of these claims for 'wet' mergers is underway, 
in this paper we would like to explore the possibility that 
such alignment (or counter-alignment) mechanism exists 
for purely gravitational interactions ('dry mergers'). In 
general we will seek to find if there is any bias in the 
individual spin distributions of black holes at close sep- 
arations (5 — 8M for starting typical full numerical evo- 
lutions) when starting evolutions with post-Newtonian 
methods at large radii with random spin orientations. 
These distributions, in turn, will then help in choosing 
configuration for full numerical simulations of close bina- 
ries. 

The paper is organized as follows. In Section [ll] we de- 
scribe the post-Newtonian formalism to analyze the in- 
spiral stage of the binary evolutions. We use the Hamilto- 
nian formulation (up to 3.5PN order) to derive the equa- 
tions of motion in the ADM-TT gauge. Conservative 
and radiative effects of the spins are included up to the 
next leading PN order. We also include a purely ana- 
lytic analysis of the projection of the quantity A along 
the orbital angular momentum L, which has a strong 
effect on the recoil velocity, to qualitatively predict a 
slight bias towards counter-alignment of these two vec- 
tors. The results of the statistics of numerical integration 
of the post-Newtonian equations of motion (EOM) fol- 
lows. We performed integrations from initial separations 
of r = 50M with 16'* spin orientation chosen at ran- 
dom and magnitudes fixed at large astrophysical values, 
i.e. Si/mf — 0.97 for different mass ratios in the range 
1/16 < q = mi/TO2 < 1- The results quantitatively con- 
firm the bias towards counter-alignment of A and total 
spin S with respect to the orbital angular momentum L. 
Section |III| deals with the merger phase, when the black 
holes are much closer to each other and in a few orbits 
will merge into a single larger one. This is the typical sce- 
nario that full numerical simulations assume. The bulk 
properties of the remnant black hole can be summarized 
in terms of empirical remnant formulae that describe its 
total mass, spin and recoil velocity. We proposed formu- 
lae for these quantities based on post-Newtonian scaling 
with amplitudes fixed by the full numerical simulations. 
With these formulae at hand, we perform statistical stud- 
ies by evaluation of these expression for random distri- 



butions of mass ratios and individual spins (magnitudes 
and spins) . Those evaluations lead to a large recoil veloc- 
ity tail in the distribution with non negligible probabil- 
ities for V > 1000 km s~*, and highly peaked about the 
direction of orbital angular momentum at merger. Like- 
wise, evaluations of the final spin formulae lead to a wide 
distribution peaked at magnitudes of Sf/M'j ss 0.73 and 
orientations peaked at an angle ~ 25° with respect to the 
orbital angular momentum. We complete the paper with 
a discussion of the astrophysical consequences of these 
results and include an appendix with the computation 
of the innermost stable circular orbit radius, energy and 
angular momentum around Kerr black holes (needed for 
the remnant formulae), with an analytic solution for the 
equatorial and polar orbits. 



II. INSPIRAL PHASE OF BHBS 
A. PN techniques 

We construct the PN equations of motion using the 
formulae provided in Refs. [54U57j . To obtain the con- 
servative part of the PN equations of motion, we use the 
following Hamiltonian, 

H = Ho ,Nowt 

+-ffsS,2PN + ^^SO,2.5PN + ^0,3PN 
+-ffSiS2,3PN + ^fsiSi(S2S2),3PN , (1) 

where Hq contains the terms associated with the orbital 
motion up to 3PN order, H^o contains the spin-orbit 
coupling terms up to 2.5PN order, and i?ss contains the 
spin-spin coupling term up to 3PN order. Note that 
Porto and Rothstein has discussed the spin-spin inter- 
action by using effective field theory techniques [58H6T] . 
These are a very powerful approach to systematically dis- 
cuss the dynamics of finite size objects. 

The equations of motion are then obtained via, 



~dr 

~dt 
dSi 

dS2 
dt 



- {X\H} = 



dH 



= {P,,H} + F, = - 



dH 



F,, 



ab2 



(2) 
(3) 
(4) 

(5) 



where {•••,•••} denotes the Poisson brackets, X"^ — x\ — x\ 
and are relative coordinates and linear momenta of the 
binary. Si and 5*2 are the spins of each body, and Fi is 
the radiation reaction force. The radiation reaction force 
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F is given by 



F = 



1 dE 



\L\ dt 



P 



15 



61 + 48 



61 + 48 
mi 



TO2 



TO2 



P-Si 



(6) 



where L = X x R = \X\, = (Muj)^/^ and w is the 
orbital frequency. We use the fohowing notation: 



M 
SM 

f] 
S 



So 




(7) 
(8) 

(9) 
(10) 

(11) 



(12) 



To calculate dE/dt, the instantaneous loss in energy, we 
use the formulae given in Refs. [52] [103] . 



1. PN prediction of distribution ofL-A 

The time derivative of the inner product L ■ A where 

L and A are the unit vector corresponding to L and A, 
respectively, is given by 



iL-A) = 



L-A 



L-A 



\L\\A\ \L\\A\ 
_L-A|L|- i-A|A| 
~ liPIAI |L||A|2 



(13) 



Here since we focus only on the dissipative effect, we ig- 
nore A and |A|' This is because there is no radiation 
reaction term in Eqs. Q and ([s]). The radiation reac- 
tion effect are introduced by the evolution equation of 
the linear momentum given in Eq. ([S]). Furthermore, we 
expect that the time evolution of the spin directions due 
to the conservative force will cancel out in a statistical 
treatment. Hence, we have 



(L-A), 



A L-A\L\- 



dis 



ILIIAI 



\mA\ 



(14) 



The dissipative effect on the angular momentum is 
given by 



^dis — 



X X F 
1 dE ^ 

T 

uj at 

-^,^4^|f61+48"^|P.5i 

15 |i|2 \V TOl.' 

+ ( 61 + 48 — ^ P • 52 P 

TO2y 



(15) 



where we have used the quasi-circular assumption and 
Eq. ([6|. 

Using this dissipation of the angular momentum, we 
obtain 



(^•A)dis = - 



1 



15 M (l + q)'^ \d2-qd1 

+q [(61 q + 48) - (61 + 48 g)] (P • ai)(P • ^2) 



-q^ (61 q + 48) (P • aif + (61 + 48 q) (P • a^Y 



(16) 



r 



in thejeading PN order calculation. Here di = Si/mj, in L^-^^ cancels out, and we have 
Q!2 ~ 82/1712 and q = mi/m2- Note that the dE/dt term 
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+q [(61g + 48)+g (61 + 48 g)](P-ai)(P -as)!, (17) 

I 

Next, we consider the time integration from i = to t = tf. 



= -^TT^ I- ^ - (61g + 48)(73.ai)2 + (61 + 48g)(P.Q;2)2 

36 (1 + g)^ |q!2 -qai\ [ 

+q [(61 g + 48) - (61 + 48 g)] (P • di){P ■ ds) 

= -^TV^ I 2-\- (61 g + 48) (^ -51)2 + (61 +48 9) (^-a^)' 

36 (1 + (j)"^ Iq^^ai + a2\ [_ 

+q [(61 g + 48) + (? (61 + 48 q)] (P • di)(P • cfs) 




3/2 



3/2 



(18) 



where we considered only the evolution of w^j, i.e., the 
inspiral, and have used the leading radiation reaction and 
the Newtonian velocity. 



dR 

Itt 



64 



(19) 



5 [l + qY 
~R ' 

In the above integration, we derived the formula as- 
suming P ■ cii ~ constant {i = 1, 2). However, since the 

spins precess, we need to check the evolution of P ■ di. 
From the evolution equations for spins in the leading PN 

order, the evolution equations for P • di are given by 



P-di] = 



P-d2 



R 



1 



P2 {l + qf 
Mv^ q 

R R^ (l + g)2 
X (X ■ d2 



2q 



(20) 



We note that in the limit g — > 0, we only need to 



consider the evolution of P • a2 in Eq- (16 1, 



(21) 



This equation means that the direction of d2 does not 
change, i.e., there is no precessio n of the spin. Therefore, 

we may replace (P • a2)^ in Eq. ( 18 1 by the one-orbit av- 
erage < [P ■ d2Y >t of (P-Q?2)^- Although the adiabatic 
evolution of < (P • c52)^ >t is present, its effect comes in 



at higher PN order in Eq. (18). In this case, it should 



be noted that we may consider a test particle orbiting 
around a Kerr black hole with the spin 5*2. According 
to |63j in the black hole perturbation approach, the par- 
ticle's angular momentum and the black hole's spin tend 
to be anti-parallel. 

On the other hand, in the case of comparable mass 
binaries, the direction of di changes on a timescale much 
shorter than the integration time. Hence, Eq. ( 18 ) is not 
expected to be accurate in g — >■ 1 limit. 



In Table |l] we show the q dependence of Eq. ( 18 1 when 
we ignore the spin precession. Here, we take the average 
with respect to the direction of two spins to represent 
the randomly oriented spins. We also present the spin 
amplitude dependence in Table In) 
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TABLE I: The q dependence in the evolution of (L ■ A)dis and 
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f T, • .S* 1 J ■ fTTiTn r 


— ^C\A4 to r — 5M Wp spt Irfi 


= Q2 = 0.97. 
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250 


1.00 


0.0000 


-0.0283 


200 


0.75 


-0.0111 


-0.0287 


0.50 


-0.0224 


-0.0310 


150 


0.25 


-0.0343 


-0.0366 




0.125 


-0.0406 


-0.0412 


100 


0.0625 


-0.0440 


-0.0441 




0.00 


-0.0475 


-0.0475 


50 



TABLE II: The amplitude dependence of the spin in the evo- 
lution of (L ■ A)dis and (L • S)dis from r = 50Af to r = 5M. 
We set q = 0.25 and = |q2| = Q. 



a (L ■ A)dis {L ■ g)dis 

0.97 -0.0343 -0.0366 

0.97/V2 -0.0242 -0.0259 

0.97/2 -0.0171 -0.0183 

0.97/4 -0.0086 -0.0092 

0.97/8 -0.0043 -0.0046 

0.97/16 -0.0021 -0.0023 



FIG. 1: The P{^i = L ■ A) distribution for g = 1 at r = 5M 
starting from a uniform distribution at r = 50A/. Here we 
plot the number of events in the given range of fj, out of 16* 
total events. 



B. Statistical Results 



For our PN evolutions with used an adaptive fourth- 
order Runge-Kutta time-integration scheme with a rela- 
tive tolerance of 10"^^. The initial data for the simula- 
tions were generated using the 3PN conservative equa- 
tions for quasi-circular orbits with orbital frequency 
Mil = 0.00275, which corresponds to an orbital radius 
of 50 ± 2M. In most cases we stopped the PN simula- 
tions at a fixed orbital radius of 5M, but also performed 
a set of simulations that terminated at r = 8M in or- 
der to see the effect of the final orbital radius on the 
distributions. To obtain the initial PN orbital param- 
eters, we used uniform distributions of q?i and a2 over 
the sphere (by choosing uniform random distributions 
in fi — cos 9 and </)) with fix amplitude a — 0.97. We 
produced 65536 random spin configurations for each fix 
mass ratio q = 1, 3/4, 1/2, 1/4, 1/8, 1/16. Each run took 
approximately 10 minutes. In addition we performed sets 
of 65536 run for q = 1/4 and a = 0.97/^2, a = 0.97/2, 
as well as a = 0.97/\/2 but terminating at r = 8M rather 
than r — 5M. We denote these three latter distributions 
in Table [m] by 0.25S1, 0.25S2, and 0.25F, respectively. 

In the following section we examine the distribution 

of the angle n — L ■ A that A — M{S2/m2 — Si /mi) 
makes with the orbital angular momentum (at r — 5M). 
At r — 50M this distribution is uniform (since Si and S2 
are chosen from a uniform distribution on the sphere) . In 
Figs. [l|6] we show histograms of the distribution of the 

angle L ■ A that A makes with the orbital angular mo- 



mentum for the given mass ratios. To analyze these data 
quantitatively, we bin the data from fi = — lto/i = 1 
with bin widths of = 0.01. We fit the resulting data 
P(yu) to a linear function P(/u) = P{0) + ^|oM for each 
mass ratio. The results are summarized in Table ITTTl and 
plots of the fits are given in Figs. 9|T4} We perform a sim- 
ilar analysis for the angle that S makes with the orbital 
angular momentum (see Fig. 17 1. We also perform a sim- 
ilar analysis, but with q fixed to g = 1/4 and ai = a2 



reduced by factors of y/2 and 2, respectively (See Figs. 7][8 
and 15p6 1, and fit the resulting slope dP/dji as a func- 
tion of a. Here the fit favors a leading-order linear de- 
pendence in a over a leading-order quadratic dependence 
(where the constant term is assumed to be zero) (See 
Fig. 18 1. If we set the constant in the fit to zero, then a 
linear dependence is dP/dn = -(0.02491 ± 0.00098)a for 

the distribution of L- A and dP/d/i = -(0.0301±0.0041)a 

for the distribution of the L ■ S. Note that the skewing of 
the distributions takes place at smaller radii, as can be 
seen by differences in the 0.25S1 and 0.25F distributions, 
which differ only in the orbital radius (5Af for 0.25S1 and 
8M for 0.25F) where the distributions are measured. In 
Table [TV] we show fits for the distributions of the angles 
Si ■ L and S2 ■ L for the same set of runs. Note that 
the distribution of • L (the smaller BH's spin) become 
essentially uniform for g < 1/4. 
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FIG. 2: The P{n = L-A) distribution for g = 3/4 at r = 5M 
starting from a uniform distribution at r = 50M. Here we 
plot the number of events in the given range of fi out of 16* 
total events. 



FIG. 5: The P{fi ^ L-A) distribution for g = 1/8 at r = 5M 
starting from a uniform distribution at r = 50M. Here we 
plot the number of events in the given range of fj, out of 16* 
total events. 





FIG. 3: The P{fi = L-A) distribution for g = 1/2 at r = 5M 
starting from a uniform distribution at r = 50M. Here we 
plot the number of events in the given range of /i out of 16* 
total events. 



FIG. 6: The P(/j = L- A) distribution for g = 1/16 at r = 5M 
starting from a uniform distribution at r = 50A/. Here we 
plot the number of events in the given range of ^ out of 16* 
total events. 





FIG. 4: The P(/x = L • A) distribution for g = 1/4 at r = 5A/ 
starting from a uniform distribution at r = 50M. Here we 
plot the number of events in the given range of ^ out of 16* 
total events. 



FIG. 7; The P{^i = L-A) distribution for g = 1/4 at r = 5M 
starting from a uniform distribution at r = 50M and cti = 
02 = 0.97/\/2. Here we plot the number of events in the given 
range of fi out of 16* total events. 
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FIG. 8: The P{fi = L • A) distribution for g = 1/4 at r = 5M 
starting from a uniform distribution at r = 50M and ai = 
02 — 0.97/2. Here we plot the number of events in the given 
range of fi out of 16* total events. 




FIG. 11: The fit to the normalized P{fi = L ■ A) distribution 
at r = 5M for q = 1/2. The data have been binned with a 
bin width of Sjj, — 0.01 
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FIG. 9: The fit to the normalized P{fi = L • A) distribution 
at r = 5M for g = 1. The data have been binned with a bin 
width of (5/i = 0.01 and normalized to a total probability of 1. 
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FIG. 12: The fit to the normalized P{fj, = L ■ A) distribution 
at r = 5M for q = 1/4. The data have been binned with a 
bin width of 5fi = 0.01 and normalized to a total probability 
of 1. 
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FIG. 10: The fit to the normalized P(^i = L • A) distribution 
at r = 5M for q — 3/4. The data have been binned with a 
bin width of 5fj, = 0.01 and normalized to a total probability 
of 1. 



FIG. 13: The fit to the normalized P(/i = L • A) distribution 
at r = 5M for q = 1/8. The data have been binned with a 
bin width of Sfi = 0.01 and normalized to a total probability 
of 1. 
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FIG. 14: The fit to the normalized P{fi = L ■ A) distribution 
at r — 5M for q = 1/16. The data have been binned with a 
bin width of 5fi = 0.01 and normalized to a total probability 
of 1. 
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FIG. 15: The fit to the normalized = L- A) distribution at 
r = 5M for (/ = 1/4 and ai = Q2 = 0.97/^2. The data have 
been binned with a bin width of S/i = 0.01 and normalized to 
a total probability of 1. 
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FIG. 16: The fit to the normalized P(/i = L • A) distribution 
at r = 5M for g = 1/4 and ai = 02 = 0.97/2. The data have 
been binned with a bin width of S/i = 0.01 and normalized to 
a total probability of 1. 



TABLE III: The distribution P{fi) of the angle = cos 61 
between A and the L at r = 5M starting from a uniform dis- 
tribution at r = 50M (top), and the similar distribution for 
the angle between S and L (bottom). The 0.25S1 configura- 
tions had ai = a2 = 0.97/\/2 and the 0.25S2 has a = 0.97/2, 
while the 0.25F configurations have a = 0.97/\/2 and provide 
the distributions at r = 8 Ad (rather than r — 5M), all others 
had Qi = Q2 = 0.97. 
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FIG. 17: The dependence of the slope in the distribution of 
the angle between A and the orbital angular momentum, as 
well as the angle between S and the orbital angular momen- 
tum as a function of mass ratio. 



An important consequence of choosing uniform distri- 
butions for the directions of Si and S2 (with magnitude 
jS'il = 0.97) is that the initial distributions for the squares 
of the magnitudes of S and A, P(5^) and P(A^), are 



uniform in the range [Q;(m| — m\)]'^ to [a{m\ + m^)]^, 
and zero outside this range (i.e. there is an equal prob- 
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TABLE IV: The distribution of the angle fi = cos 61 

between 5*1 and the L at r — 5AI starting from a uniform dis- 
tribution at r = 50M (top), and the similar distribution for 
the angle between S2 and L (bottom). The 0.25S1 configura- 
tions had Qfi = 0:2 = O.97/V2 and the 0.25S2 has a = 0.97/2, 
while the 0.25F configurations have a — 0.97/\/2 and provide 
the distributions at r = 8M (rather than r = 5M), all others 
had ai = 02 = 0.97. Note that the distribution of angles for 
the smaller component Si becomes uniform as g — > 0. 
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ability of finding any given value of S*^ or in this 
range). However the distributions P(A) and P{S) there- 
fore contain an additional linear factor in A and S (i.e. 
P{x) = 2xP{x^) for any variable cc), respectively. One 
immediate consequence is that the distributions P(A) 
and P{S) are both maximized for the largest allowed 
values of S and A. Given the observation that large A 
in the orbital plane [J leads to very large recoils, this 
bias, if present in nature, would favor observations of 
large recoils. See Sec. |IV| for further analysis of the recoil 
distribution. 

Schnittman in Ref. ^645 has studied the evolution of 
spins in binary systems using orbit-averaged PN equa- 
tions of motion what allowed longer term evolutions 
(from separations up to lOOOM). The results indicate 
strong correlations of the late angle among spins when 
one starts fixing the initial direction of the spin of the 
primary object and choose the secondary's spin direction 
at random (See Figs. 6 and 7 in [H].) Bogdanovic et 
al revisit this scenario in Ref. [511 and find that if one 
is allowed to choose initial random distributions for both 
spins the resulting evolution leads to close to isotropic 
distributions of the late directions of the spins (See theirs 
Fig. 1). In our paper wc find an small but statistically 
significant bias towards counteralignment of the spins 
with the orbital angular momentum (See Figs. 9p6 ) 

More recently, Herrmann et al. presented numerical 
studies of the PN equations on GPUs [5S]. They used 
the evolution equation for the orbital frequency assum- 
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FIG. 18: The dependence of the slope in the distribution of 
the angle between A and the orbital angular momentum L for 
q = 1/4 as a function of |(5i| = |a2| = a, as well as the angle 
between S and the orbital angular momentum as a function 
of a. In all fits the constant term is taken to be zero. The 
data here favor a linear dependence in a. 



ing quasi-circular orbits. This equation is coupled with 
the spin and angular momentum precession equations, 
which include the leading order spin-orbit and spin-spin 
couplings. On the other hand, in our calculation, the 
PN equations of motion are derived from the Hamilto- 
nian and include radiation reaction effects. These have 
higher PN order spin-orbit and spin-spin coupling terms. 
Furthermore, the second term of the right hand side of 
Eq. (|6| has a significant effect in the PN evolutions. Al- 
though the evolution of L in |65j is determined only by 
the conservative dynamics, we have also considered the 
dissipative effect due to the radiation reaction. We find in 
the PN prediction that this dissipative effect creates the 
statistically significant counter-alignment of the spins. 



III. MERGER PHASE OF BHBS 

Unlike in the earlier inspiral phase, during the plunge 
and merger the PN equations of motion do not provide 
a quantitatively accurate description of the merger dy- 
namics, and therefore do not provide robust estimates of 
the final remnant mass, spin, and recoil. However, analy- 
sis of the recoil in particular shows that PN analysis can 
be used to derive heuristic formulae (based on how PN 
predictions scale with spins and masses) that give quan- 
titatively correct predictions ,7, ,24. 166) and incorporate 
the symmetries of the problem. We will use this model- 
ing in the case of the total radiated energy and angular 
momentum. In particular we will supplement the inspiral 
losses, modeled by the energy and angular momentum of 
the ISCO in the particle limit (extended to the compa- 
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rable mass regime) with the subsequent plunge using the 
PN dependence on the BHs parameters (and fitting the 
amphtudes as in the recoil velocities case). 



A. Recoil velocities 

In order to quantify and model the nonleading cor- 
rections, we augment our original empirical formula 
with new subleading terms that are higher order in 
the mass ratio and include a new term linear in the 
total spin, motivated by higher order post-Newtonian 
computations [67' , and introduce additional parameters 
Bh,Bk,Hs,Ks and Gi, 



K-ccoii(g, a) 



Vm ei + v± (cos ^ ei -I- sin C 62) + v\\ n\\ , 



A 



H 



(1 + 9) 



[l + Brj] 



(1 + 9) 

(1-9) 



(1 + 9) 



(1 + Bh v) {0^2 - lal) 
2 (a| + q^al) , 



= K- 



(1 + Bk v) - ic'i 



(1 + 9) 

X cos(9a — Bo) 



X cos(85 — 61) 



(22) 



where rj = q/{l + q)^ , with q = mi/m2 the mass ratio of 
the smaller to larger mass hole, di = Si/mf, the index 
_L and || refer to perpendicular and parallel to the or- 
bital angular momentum respectively, ei , 62 are orthogo- 
nal unit vectors in the orbital plane, and ^ measures the 
angle between the unequal mass and spin contribution 
to the recoil velocity in the orbital plane. The constants 
Hs and Ks can be determined from newly available runs. 
The angle 8 is defined as the angle between the in-plane 
component of A = M{S2/m2 — Si/mi) or S = Si + S2 
and the infall direction at merger. Phases 9o and 81 
depend on the initial separation of the holes for quasicir- 
cular orbits. 

A crucial observation is that the dominant contribu- 
tion to the recoil is generated near the time of formation 
of the common horizon of the merging black holes (See, 
for instance Fig. 6 in [SS])- The formula above (22) 



describing the recoil applies at this moment (more pre- 
cisely, the coefficients correspond to an averaging of the 
PN expressions during the plunge phase), and has proven 
to represent the distribution of velocities with sufficient 
accuracy for astrophysical applications. The total recoil 
velocity also acquires a correction [52] for small eccen- 
tricities, e, of the form Ve = Kccoii (1 + e), and if one 



allows for relativistic close hyperbolic encounters, then 
recoils up to 10000 km s~^ are possible [70]. Although 
we expect the orbits will circularize well before merger. 

The most recent estimates for the above parameters 
can be found in |50^ and references therein. The current 
best estimates are: A = 1.2 x 10^ km s^^, B = —0.93, 
H = (6.9 ± 0.5) X 10^ km s"\ K = (6.0 ± 0.1) x 
10* km s~^, and ^ ~ 145°. Note that we can use the data 
from [50] to obtain K = (6.072 ± 0.065) x 10'' km s~\ 
if we assume that Bk and Ks are negligible. On the 
other hand, by fitting the data to K and Bk simul- 
taneously, we obtain K = (5.24 ± 0.29) x 10* km s^^ 
and Bk = 0.74 ± 0.29. At first glance these two re- 
sults look quite different. However, in both cases the 
actual resulting empirical formula predict the same re- 
coil velocities within 8% over the range 1/10 < g < 1 
(the data from [50] covered the range 1/8 < q < 1). 
Finally, if we fit the data to find K and Ks simultane- 
ously we obtain K = (6.20 ± 0.12) x 10* km s^^ and 
Ks ~ —0.056 ± 0.041, where we made the additional as- 
sumption that since S* = A for these runs, that 80 = 81. 
An attempt to fit all three parameters produces inaccu- 
rate fitting parameters because the degrees of freedom in 
the fit and the limited number of available runs. 



Equation (22 1 for the recoil contains all the expected 



linear terms in the spin, and include ten fitting parame- 
ters. Based on the works [67j[7T] one could add quadratic 
terms, but they are complicated expressions with more 
fitting parameters that we will not include here. 



B. Remnant Mass 

Motivated by the success of the empirical formula for 
the recoil, we propose a new empirical formula for the 
total radiated energy based on the post-Newtonian equa- 
tions that describe the instantaneous radiated energy 
(See Eqs. (3.25) in Ref. [55], and for the quadratic terms 
in the spin see Ref. [57], Eq. (5.6)). For example, the 
spin-spin contribution to the radiated energy has compo- 
nents quadratic in A that have the form. 



E 



ss AA'^ + B{h-Af + C{v-Af+D{n-A){v-A) 

= A{A]_+Af) + A]_{I3coa^e + Ccoaesme + D) 

= AAl + bAl(cos^{0-9o) + c) 

= AA?:+bA']_{cos2{e -eo)+£). (23) 



A similar expansion can be derived for the terms 
quadratic in 5^o = 25 -|- {SAI/M)A. In addition to the 
terms arising from the instantaneous radiated energy, 
that allow for twelve fitting parameters, we also included 
terms associated with the secular loss of energy in the 
inspiral period from essentially infinite separation down 
to the plunge. In order to model this contribution we 
adopted the form of the the 2PN binding energy, with 
coefficients chosen to reproduce the particle limit at the 
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ISCO[72] 



Eisco = (l-V8/3) + ^ 



324%/2 



where we considered terms up to quadratic order in the 
spin. 

If we take into account the 7/^ effects from self force 
calculations [73J and 2PN effects of the spins (See [i^ . 
Eq. (4.6)), we obtain: 



E 



isco 



(1 - \/8/3) + 0.103803?? 
1 



36^/3(l + g)2 
5 



q{l + 2q)a\ + {2 + q)a}^ 

_ (So — 3(q;!)^ 

32472(1 + 9)2 ' ^ 2^ 

~2q(di ■ d2 — 3af q;|) + q^iaf — 3(af )^) 

+0{a^) (24) 

This expression represents a quadratic expansion in the 
spin-dependence, hence we expect to produce reliable re- 
sults for intrinsic spin magnitudes < 0.8. The exact 
expression for all values of spins, including maximally 
rotating, are complicated and are given in the appendix 
Thus our parametrization of the energy loss is given 

SM/M = Eisco + -^2??^ + £^3??^ 

2 f 

2 



+ 



Es{a!^ +q al) 



(1 + 9)' ^ 

+Ea (1 - q) (a| - qal) + Ea |a2 + ^^iP 
+Eb \ai+qai\^ (cos2(e+ - 62) + Ec) 
+Ed \d2 - qaip 



^Ee K 



qai\' 



W(e_ -e 




where M = mi + m2 and Q± are the angles that 
A± = M (Si /mi ±62/7712) make with the radial direction 
during the final plunge and merger (for comparable-mass 
BHs, a sizable fraction of the radiation is emitted during 
this final plunge, see for instance Fig. 6 in Ref. |68)). 
Phases 82,3 are parameters that give the angle of maxi- 
mum radiation for these terms, and depend on the initial 
separation and parameters of the binary at the beginning 
of the numerical simulation. 

According to the FN theory [Ti] , the leading correction 
to the radiated energy for small eccentricities has the 
form E = £'c(l + 157/24e), where Ec is the radiated 
energy in the circular case, which should in principle be 
added to the above formula. However, it is expected that 
the orbits will be quite circularized by the time of merger. 

To determine the fitting parameters in formula ( 25 1 we 



for energy already lost by the system in reaching the ini- 
tial separation of the simulation. That is, some authors 
choose to normalize their data such that the sum of the 
horizons masses is 1, which approximates the situation 
where the energy lost during the prior inspiral is taken 
into account, while other normalize the initial data to 
unit ADM mass. In these latter cases, we add the 3FN 
binding energy of the initial configuration to the calcu- 
lated radiated energy to obtain an approximation for the 
total energy radiated by the system in question from in- 
finite separation. 

For the non-spinning case we fit the data found in 
Refs. [751 ES]- Here we fit i^Rad versus 77, where £^R,ad 
is the total radiated energy for a given configuration mi- 
nus the binding energy of the initial configuration (where 
the binding energy is negative) . We calculate the binding 
energy using the 3PN accurate expressions given in [77] . 
A fit of the resulting data gives E2 = 0.341 ± 0.014 and 
E3 = 0.522 ± 0.062. 

For the spinning cases, when spins are aligned with the 
orbital angular momentum, fits for final remnant mass 
from Ref. 78, 79 yield Es = 0.673±0.035, E^ = -0.36± 
0.37, Ea = -0.014 ± 0.021, and Ed = 0.26 ± 0.44. The 
source of these large errors is the difference in correcting 
for the normalization of the results in the papers. 

Finally, fits from the final remnant masses from 
Ref. [24| yields Ee = 0.09594 ± 0.00045 and fits from 
the equal-mass configurations in Ref. [5D] yield Eb ~ 
0.045 ±0.010. 



C. Remnant Spin 

In an analogous way, we propose an empirical formula 
for the final remnant spin (note that the total radiated 
angular momentum, unlike the total radiated energy, is 
not finite for an inspiral from infinite initial separation) 
based on the post-Newtonian equations that describe the 
radiated angular momentum (See Eqs. (3.28) in [IS] ) 
and the angular momentum of a circular binary at close 
separations (4.7), 



ttfinal — 



(1 - SM/M)-^ {vJisco + + hri^) n\\ 



2 

+ (T+^(['^^^"2+9'«!) 

+Jb (1 - q) (a| -qa{) 

+ {l-q)\ai -qdi\ 

X 7 Ja cos[2(eA - 64)] + Jma 7J_L 
+ \a^+q^ai\ 

X ^Js cos[2(es - 65)] + Jms }• 



(26) 



need to correct some of the numerical data to account 



where we have expanded the triple cross products in Eq. 
(3.28c) in [46J and used the last form of the parametriza- 
tion in Eq. 



13 



Note that, even at linear order, there are important 
contributions of generic spinning black holes producing 
radiation in directions off the orbital axis that do not 
vanish for equal masses nor vanishing total spin. The 
above formula can be augmented by quadratic-in-the- 
spins terms [67l [71] of a form similar to the terms added 
to the radiated energy formula (25). However, those 



terms are more complicated and involve many more fit- 
ting constants in addition to the ten for the linear depen- 
dence. In addition, the linear approximations seem to 
have smaller quadratic corrections for the radiated angu- 
lar momentum than the radiated energy (See for instance 
Fig. 21 of Ref. ^.) 

Looking at the spin expansion of the orbital angular 
momentum of a particle at the ISCO [72] . 



from the data in [75] and [73], we find J a = -2.97 ± 0.26 
and Jb — —1.73 ± 0.80. From the combined fit we find 
that 2.42% < 5M/M < 9.45% and 0.34 < a < 0.92 for 
the equal-mass, aligned spin scenario, in the region where 
the fit is valid (|a| < 0.9). 



Note that formulae ([24j) and ( 28 1 generalize the particle 
limit ISCO results to take into account both the mass 
ratio dependence, q = mi/m2, and the spin of the holes. 
Si and S2, in a symmetric way that is accurate up to 
quadratic order in those parameters. This allows us to 
have an accurate description both when the binaries have 
relatively small mass ratios, i.e. q ^ 1/10, and in the 
comparable mass regime (where the radiative terms in 
(251 and (26) dominate and the ISCO is ill defined). For 



the full expressions see Appendix A. 



Jisco = 2\/3nll - 



9\/2 



a2 



9^3 



-,2 

"2 



a2 



r;(l + q)2 



(27) 



and incorporating the ry^ effects from self force calcula- 
tions [73] and the 2PN effects of the spins (see Eq(4.7) in 
[46]). fitted to reproduce the particle limit one obtains, 

Jisco = < 2^3 - 1.525586277 



q{7 + 8q)al + (8 + 7(7)a| 
«2-3(a|)2 



9^/2(1 




2 




9^3(1 




2g(ai • 


0L2 - 


1 




9^/2(1 


+ 9)2 


1 ((32 -i 


-q^oLx 


V (1- 





^2/^2 



[q{l + 4(7)ai + (4 + q)d2] 
-+0(a3). 



(28) 



This expression represents a quadratic expansion in the 
spin-dependence, hence we expect to produce reliable re- 
sults for intrinsic spin magnitudes ai < 0.8. The exact 
expression for all values of spins, including maximally 
rotating, are complicated and are given in the appendix 

According to the PN theory [72] , the leading correction 
to the radiated angular momentum for small eccentrici- 
ties has the form J = Jc(l -|-23/8e). This correction can 
be added to the results, but we expect very low eccen- 
tricities by the time the BHs merge (see [101 IHI] for the 
effects of eccentricity on the remnant). 

For the non-spinning case we fit the data found in 
Refs. [751175]. We find J2 = -2.81 ± 0.11 and J3 = 
1.69 ± 0.51. Fits for final remnant spin from Refs. [78] 
yield Ja = -1.971 ± 0.018, and Jb = -3.611 ± 0.042. 
The rest of the fitting constants are currently hard to 
determine with precision. If we attempt to fit Ja and Jb 



IV. MERGER STATISTICS 

In order to predict the distribution of the remnant spin 
and recoil, we consider a normal random distribution of 
spin directions and magnitudes and mass ratios for quasi- 
circular black-hole binaries and compute the probability 
density distribution of spin magnitudes of the final rem- 
nant (see also [20, 82J). Although our initial PN simula- 
tions showed a slight bias towards counter-alignment of 
the spins to the angular momentum, the effect is small 
and will be neglected here. The results of ten million 
of such simulations are displayed in figure 19 The fi- 
nal results are insensitive to the initial distribution and 
quickly converge, in a few generations of mergers, to the 
displayed curve, which consequently represents a univer- 
sal distribution of the intrinsic spin magnitudes [with a 
maximum near 0.73 and mean in the range (0.53, 0.83)] 
of remnant BHs (when spin-up effects due accretion are 
not taken into account). We fit the distribution of the 
final spins to the Kumaraswamy functional form |83j 
fix; a, b) = a5a;°-i(l - a;°)''"\ and find a = 6.58 ± 0.08, 
& = 7.14 ±0.19. 

Figure [20] show the probability distribution of the mag- 
nitude of the final remnant hole's spin for different ranges 
of the mass ratio q. We observe that the distribution 
becomes more and more peaked around the value of 
a ~ 0.75 as we go to comparable masses. The distri- 
bution gets wider for small mass ratios, this is expected 
because in the g — )■ limit one get essentially the origi- 
nal spin distribution of the massive black hole that was 
taken randomly. Compare this plot with a similar one in 
Ref. [5D], Fig. 1, for dry mergers. 

Figure [22] displays the angular dependence of the prob- 
ability distribution of the final spin with respect to the 
original orbital angular momentum at far separations for 
different ranges of mass ratios. The distribution for com- 
parable masses is peaked at angles close to the orbital 
angular momentum since the spin contributions tend to 
cancel leaving most of the contribution to the final angu- 
lar momentum to the orbital component. We also observe 
that as we go to smaller mass ratios the distribution be- 



14 



P(a) 




FIG. 19: The remnant BH spin magnitude distribution after 
several generation of mergers, starting from a uniform initial 
distribution. 



P(ff) 




FIG. 20: The remnant BH spin magnitude distribution after 
several generation of mergers (starting from a uniform ini- 
tial distribution) for mass ratios in the ranges < 5 < 0.1 
(dashed), 0.1 < g < 0.2, 0.4 < g < 0.5, 0.6 < g < 0.7, and 
0.9 < q< 1.0. The distribution becomes more sharply peaked 
around larger values of a as the mass ratio increases. Note the 
for the smallest range, the spin is more highly peaked than 
for 0.1 < g < 0.2. 
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FIG. 21: The remnant spin direction distribution after several 
generation of mergers, starting from a uniform initial distri- 
bution. 
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FIG. 22: The remnant spin direction distribution after several 
generation of mergers, starting from a uniform initial distribu- 
tion, for mass ratios in the ranges < g < 0.1, 0.2 < q < 0.3, 
and 0.9 < q < 1. The closer to equal masses, the more 
highly peaked the distribution. For mass ratios in the range 
0.9 < q < 1, the distribution is peaked at 6 ~ 25°. For the 
smaller mass ratios, the distribution approaches sin 61 (i.e. the 
uniform spin direction distribution). 



comes wider since the larger black hole contributes ran- 
domly to the final total angular momentum. We also 
observe the vanishing probability of having exact align- 
ment of the final spin with the initial angular momentum. 
This is because exact alignment is a set of measure zero 
on the initial random distribution of the spins. 

We then use these distributions for the spin- 
magnitudes, while assuming uniform distributions in an- 
gle and mass ratio, to predict the distribution of recoil 
velocities. In Figs. [24] and [27] we plot the distribution 
of recoil magnitudes and directions, respectively. From 
the plots we see that distribution of recoils in angle is 
strongly peaked toward alignment and counter-alignment 
with the orbital angular momentum, which also gives the 
largest recoil magnitudes (this is a consequence of the fact 



that the out-of-plane recoil is generally an order magni- 
tude larger than the in-plane recoil, as seen in Eq. ( [22] ) 
and any small component of the spins along the orbital 
plane lead to those large off-plane recoils.) In Table V] 
we show probabilities for producing various large recoils. 



In Fig. 23 we plot the distribution of recoil velocities as 
a function of the angle that the recoil makes with the or- 
bital angular momentum. Here we find that most recoils 
are aligned (or counter aligned) with the orbital angular 
momentum (the distributions of recoil velocities P{v) at 
an angle 9 with respect to L is identical to the distribu- 
tion P{v) at an angle tt — 9). This strong dependence of 
the magnitude of the recoil with the angle off the orbital 
plane was to be expected given the large anisotropy found 
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TABLE V: The probability to obtain large recoil velocities, 
and large recoil velocities along the line of sight. 
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FIG. 23: Velocity distribution as a function of the angle the 
recoil makes with the orbital angular momentum (the distri- 
bution is symmetric with respect to d ^ tt — 9). The plot 
shows P{v) versus v for recoils with angles in the ranges 
(0,10°) (top), (10°, 20°) (2nd), (20°, 30°) (3rd), (80°, 90°) 
(bottom). The magnitude of P{v) is the total number of 
simulations with recoil-velocity directions between the given 
angles and magnitude within the range v ± 15 km s-^ The 
maxima occur at v(km s~^) ~ 600, 300, 200, 100, respectively. 
These distributions were obtained starting from binaries with 
a uniform distribution in mass ratio and distribution in spin- 



magnitude (random directions) given in Fig. 19 



in the empirical recoil formula ( 22 ) where the off-orbital 
plane velocities are an order of magnitude larger than the 
in-plane ones (i.e. the values of the fitted constants H 
versus K.) 

Figure |26] displays the detailed dependence of the dis- 
tribution of the magnitude of the recoil velocities with 
the mass ratio q. Here we observe that for comparable 
masses we obtain a long tail of large velocity probability 
which recedes towards smaller velocities as we reduce the 
mass ratio. This has to do with the suppressing factor 
This behavior like 



1]^ in Eqs 



(22) 



q for small mass 
ratios also explains why the probability shows a peak 
around v = 0. If we consider the probability density func- 
tion P{v) then dP = P{v)dv = P{v).{dv / dq) dq. Hence 
P{v) = {dP / dq) / {dv / dq) , and since we have chosen a 
white distribution of q we have dP/dq = 1 we obtain 
P{v) — l/{dv/dq) ^ 1/q, where we have taken the lead- 
ing dependence from Eq. ( |22[ ) w ~ q^. This explains the 
sudden growth of the probability near v = when we 
consider the velocity distribution in the mass ratio range 
< g < 0.1. The wide distributions at intermediate 
mass ratios has to do with the additional dependence on 
the direction of the spins of the holes. This distribution 
drops to near zero again near the maximum recoil veloc- 
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FIG. 24: The recoil velocity magnitude distribution for a 
uniform distribution in mass ratio and spin-magnitude dis- 
tribution in Fig 19 (with uniform spin direction). Here 
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FIG. 25: The distribution of the recoil velocity along an 
observers line of sight for a uniform distribution in mass 
ratio and spin-magnitude distribution in Fig 19 (with uni- 
form spin direction). Here < v >— 315 km s~^ and 
^/<v'2 > - <v>'2 = 358 km s'^ . 



ities V > 3000 km s~ since this configurations require 
not only comparable masses, but also counteralignment 
of near maximal spins that contribute with a set of mea- 
sure zero to the total probability density. 



V. DISCUSSION 

In this paper we studied the 'dry' (i.e. gravitational 
radiation driven) inspiral of black-hole binaries. We per- 
formed the evolutions from separations of the order of 
50M, using up to 3.5PN accurate expressions that rep- 
resents an excellent approximation for this regime (simi- 
lar simulations with full numerical relativity could prove 
practically impossible with current technologies). The 
statistical results show a small bias towards counter- 
alignment of the vectors A and S with respect to the 
orbital angular momentum L just prior to merger. This 
effect essentially takes place at close separations and can 
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FIG. 26: The recoil velocity magnitude distribution for a uni- 
form distribution in mass ratio and spin-magnitude distribu- 
tion in Fig 19 (with uniform spin direction). The plot shows 
the recoil velocity distribution for mass ratios in the range 
< q < 0.1, 0.1 < g < 0.2, 0.3 < g < 0.4, and 0.9 < g < 1. 
The distributions become successively broader for larger val- 
ues of g (i.e. similar masses). 
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FIG. 27: The recoil velocity direction distribution. The angle 
9 is in radians, and is defined as the angle between the recoil 
vector and the orbital angular momentum vector. 



be studied analytically at low post-Newtonian orders. 
The antialignment effect is associated with the late-time 
precession of the orbital plane due to radiation reaction. 
This effect for 'dry' mergers seems to oppose the align- 
ment mechanism observed in 'wet' mergers {51L I52j . 

After the initial inspiral regime, we studied the merger 
of black-hole binaries using full numerical simulations. 
Here we provided a framework to describe the bulk prop- 
erties of the remnant of a BHB merger based on PN scal- 
ing with free parameters fixed by fitting the results of full 
numerical simulations. We have shown how to determine 



the mass loss in each BHB encounter ( 25 ) and the spin of 



the remnant BH ( 26 ) with fitting constants set to match 



currently available runs. Using the same methods, one 
can improve the fitting parameters in the above formulae 
as the results of new runs are made available, thus provid- 
ing a standard method to incorporate all full numerical 



results into astrophysical modelings. Precessing, highly 
spinning binaries, including those with small mass ratios, 
simulations can provide information to also fix the non- 
leading terms in the fitting formulae. However, we expect 
that these extra terms are relatively small and will not 
have a significant impact on the results presented here. 

The new formulae are physically motivated, as they are 
derived using the post-Newtonian behavior, and natu- 
rally incorporate the correct mass ratio q dependence and 
physical symmetries, as well as allow for the radiation of 
angular momentum in the orbital plane. These formulae 
model the final plunge of comparable masses black holes 
in an impulsive approximation are supplemented by the 
slow inspiral losses that precede this regime by adding 
the ISCO energy and angular momentum in the particle 
limit, generalized to symmetric dependence on the mass 
ratio and spins (Eqs. (24) and (28)). 

We also extended the successful recoil formula by 
adding nonleading terms to include all the linear depen- 
dence with the spins, as well as higher mass ratio powers 
in Eq. ( 22 ) . Unlike in the formula for the remnant recoil 



case, the energy and angular momentum lost by the bi- 
nary during the inspiral phase is a non-trivial fraction of 
the total radiated energy and angular momentum (and, 
in fact, is the dominant contribution in the small mass 
ratio limit). We thus included both the instantaneous 
radiative terms for the plunge phase as well as the bind- 
ing energy and angular momentum at the ISCO into our 
empirical formulae (25 1 and (26). 



Using the fitted coefficients in the above formulae, we 
find that for equal-mass, non-spinning binaries, the net 
energy radiated is 5% of the total mass and the final spin 
is a ~ 0.69, both in good agreement with the most ac- 
curate full numerical runs [51]. For maximally spinning 
BHBs with spin aligned and counter-aligned we estimate 
that quadratic corrections lead to radiated energies be- 
tween 10% and 3% respectively. As for the magnitude of 
the remnant spin, the linear estimates are between 0.97 
and 0.41 respectively, with quadratic corrections slightly 
reducing those values. These results show that the cosmic 
censorship hypothesis is obeyed (i.e. no naked singulari- 
ties are formed) and are in good agreement with earlier 
estimates [5]. 



The set of formulae ( 25 ) and ( 26 ) with the fitting con- 
stants determined as in the Sec. IIIII can be used to de- 
scribe the final stage of binary black holes mergers in the- 
oretical, N-body, statistical studies in astrophysics and 
cosmology (STJ |551 155H55] by choosing a distribution of 
the initial intrinsic physical parameters of the binaries 
{q, Si, S2) and mapping them to the final distribution of 
recoil velocities, spins and masses after the mergers. Here 
we performed initial studies and have found that: i) The 
merged black holes have a considerable probability (23%) 
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to reach recoil velocities above 1000 km s~ (See Fig, 
and Table |v| and the distribution is highly peaked along 
the orbital angular momentum (See Fig. 27). ii) The 



direction of the spin of the final merged black holes is 
strongly peaked at an angle of w 25° with respect to the 



17 



TABLE VI: The probability to obtain large recoil velocities, 
and large recoil velocities along the line of sight if the distribu- 
tion in mass ratios is uniform in log^Q q for —2 < logjQ q < 0. 
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0.003% 



orbital angular momentum pre-merger (see Fig. 21 ), and 
the spin magnitude is strongly peaked at S'//MJw 0.73 
(see Fig. 19 1. Higher spins are likely if we include the 
effects of accretion. This information can be useful in 
modeling the observational effects of supermassive black 
holes kicked out of their host galaxies [96]. For example, 
as a first approximation, one may assume that the inner 
accretion disk is associated with the orbital plane of the 
merging binary, while the direction of the final spin is as- 
sociated with the current direction of the radio-jet, and 
finally that the preferred direction of the kick is along 
the orbital angular momentum. We can then to recon- 
struct 3D recoil velocities out of the observer (redshift 
velocities) information. Also, when modeling the effects 
of kicks on accretion disks surroundings the merged bi- 
nary, one should take into account that the most likely 
recoil velocity depends on the angle with respect to the 
binary's orbital plane (See Fig. 



23). 



In order to take into account that one of the main 
methods to search for recoiling black holes is to look for 
large differential redshift, typically between narrow band 
emission lines coming from the host galaxy and broad 
emission lines from the portion of the accretion disk that 
the recoiling black hole carries with it |551 - HT] . we have 
computed the effect of projecting the computed recoil 
velocity from the merger of two black holes along the 
line of sight of an observer on earth. The results are 
plotted in Fig. [25] and in Table [V] 

Finally, we note that the recoil distributions are sensi- 
tive to the assumed distribution of mass ratios. Ideally 
one should use a distribution consistent with the true dis- 
tribution of mass ratios of merging galaxies. In |38j this 
distribution is derived analytically from merger scenar- 
ios in a cosmological model. They find that the distri- 
bution is nearly fiat in logj^p 1 from q = \/ 100 to g = 1 
(see Ref. [38], Fig 1.). When we use a distribution uni- 
form in logj^Q q from —2 < logj^Q q < here, as suggested 
in _3^ , the expected recoil velocity distribution is skewed 
towards much lower velocities. In Fig. 28 we plot the re- 
coil magnitude distribution for this distribution, while 
in Table |VI| we give the probabilities for obtaining large 
recoil velocities with this mass ratio distribution. Note 
that this distribution is, in fact, strongly skewed towards 
lower mass ratios when compared to the distribution uni- 
form in q. Consequently we see much lower probabilities 
for large recoils. 



FIG. 28: The distribution of recoil velocity magnitudes and 
directions (with respect to the orbital plane) assuming a dis- 
tribution in mass ratios uniform in log^g 9 in the interval 
-2<logiog<0. 
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Appendix A: Innermost stable circular orbit of 
"Kerr" geodesies 

In this appendix we provide the necessary formulae to 
compute the Ejsco and Jisco denoted in the empirical 
equations for the black hole remnant mass and spin, i.e. 
Eqs. (25) and (26). Note that in the main text we have 



given explicitly those functions up to quadratic terms in 
the spin. We also provide explicit analytic expressions of 
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the ISCO radius for equatorial and polar orbits (See Eq. 
( |A2T| ).) 

When we treat the spin effects on the equations of 
motion in the effective one body (EOB) approach, it is 
usef ul to define two combinations of the spins |97j . Sq in 
Eq. ( 12 ) and S, 



S = 

_ 1712 

mi 



SM 
Id 



Si 



A 

mi 
m2 



^2, 



(Al) 



This is because we can rewrite S by using the nondimen- 
sional spins, c?i and q?2j 



S 



mim2 (tti + 02) , 



(A2) 



and then when we consider 77 = 0, mi or TO2 is zero. We 

note that this definition of S is same as cr in [97] except 
the numerical coefficient. 

Based on the equations given in [55] , we focus only on 
So to derive the innermost stable circular orbit (ISCO) 
of the Kerr spacetime with a spin a = Sq/M. Here, 
we assume that the direction of Sq along the 2;-axis and 



we use the Boyer-Lindquist coordinates. In practice, we 
should use and write "spherical" orbits due to the spins 
of a binary. But here, we call "circular" orbits. 

When we consider the geodesic motion, z"(t) = 
{tz{T),rz{T),9z{T),(l)z{T)}, where r is the proper time 
along the orbit, there are three constants of motion as 
follows. 



E 



(A3) 



where = dz^/dr, and the two Killing vectors are 
^^^^ = {dtY and ^1*^^ = {d^y\ And also, we define 

the Killing tensor, i^T^^ = 2SZ(^nj,) + r^gf^i,, where = 
(r2 + a^. A, 0, a) /A and n'' = (r^ + a^, -A, 0, a) /{2T) 
are two radial null vectors. Here, A = — 2Mr + and 
S — + cos^ 6. The Killing tensor satisfies the equa- 
tion K(^^^.p) — 0. We also define another notation for the 
Carter constant, C — Q — {aE — L^)^. For equatorial 
plane orbits, C defined by this vanishes. 

The energy per unit mass for circular orbit with the 
radius tbl is derived as 



E = 



2 Ma (rBL^-2MrBL 



»"BL 



^BL 



yaHP + 2ya^rBL -QyrBifa^M^ -Aya^MrBif -Qr^iI'My + r^i^^y 



-dM^BL^y ~4aH4rBL^ -QtbiI'M + tbl' + OM^rBL^) {a^V + 2rBL^ya^ + vbl^ + VBL^y)) 



hr-BL (2rBL y + rBL V 



9„,2 



23 rBL^Mya^ - 12 vbl'^ Ma'^y + 33 M^rBL^ya^ + 14 M'^rBL^a^y - 20 VBL^y^Ma^ - 18 VBL^y'^Ma 



2 rBL^a'^y + 6 rBh^y^a"^ + 4 vbl^ ya^ + 4 vbl^ y^a^ + tbl a"y 

f2 



8. .2 



+28 M^BL^y^a' + 8 M^BL^y'a^ - 4 rBL^^'a^M 



Af^BL^ya' - 4AfVBL'a^2; -iaHlrBt'y 



+aSMVBL y - 8 M^TBi^y'^a^ + 4 M^rBi.^y'^a^ ~ 4 MVbl - 3 tbl^o^ Af - 14 VBiI^My 

+5 A^Vbl^o' + 32 M^BLV - 7 tbl'^i/M + 16 A/Vbl'^' + 4 rBL^y'a' - 24 M^rBL^'y - 12 M\bl%^ 

+a^y^M ~ Ttbl'M + IGM^bl - 12 Af Vbl"^) / ((2ya4AfrBL + ^oVbl' + ya^M^ + 2yaVBL^ 

-6 yrBL^a'M^ - 4 ya^Af ^g^' - 6 tbl' A^y + rBL^y + 9 Af^rBL^y - 4 a^Af ^bl' - 6 tbl' Af 

-,1/2 

+rBL' + 9 A/Vbl^) (a*y + 2 rBL^ya' + tbl^ + ^BL^y)) • (A4) 



Here, y is introduced as a dimensionlcss inclination pa- This is related to an inclination angle as 
rameter defined by ^ 

cos iBL = , — — r 

vy + 1 

= , , (A6) 

C / A c\ the inclination angle gives the exact inclination in 

^-72- ^^"^ the case of the Newtonian orbit. 
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It should be noted that although we can define the 
circular orbit as a orbit with a constant radius in the 
Boyer-Lindquist coordinates, the circular orbit in another 
coordinates has a time dependent radius [BT . The detail 
analysis of the gauge transformation have been discussed 
in [M]- 

And then, the angular momentum per unit mass along 
the z-axis is calculated by 



Here and hereafter, we focus only on the case that > 0. 



Li 



The ISCO radius in the Kerr spacetime is obtained by 
(A7) solving the following equation with respect to tq. 



= (-6ro^Af + 4a2Mro^-6a^Mro + aS + 3aVo'^ + 3aVo2+ro^)%2 

-2 ro'' (3 - 12 tq a^M + 8 r^^a^ + 28 a^Af - 60 ra^Ma^ + 6 a^ro^ + 24 A-f Vo^'a^ 

+28 rii'a^M - 36 M^n^ + 12 tq^M - r^^) y 
+ro^ (9 - 28 o^Mtq - 6 ro^a^ + 36 M^ro^ - 12 Mtq^ + tq^) . (A8) 



In the following, we discuss the equatorial and polar 
orbits analytically. In general inclined orbit cases, we 
need to solve Eq. ( A8 ) numerically. 



1. Equatorial circular orbit 

For the equatorial orbit, we may consider y = in 
Eq. jASl). 



= (9a4 - 28a^Mro - 6ro2a2 +36M2ro2 
-12Afro3 + ro4) 

^3 - ro^ + 8 \/A7^a + 6 Afro , 



X (^3a2 - ro^ - 8VMv?^a + 6Afro) . (A9) 

This case has been discussed analytically in J 00 , which 
we summarize below. It should be noted that we consider 
the parameter range —M < a < M and only the orbits 
with > in our treatment. 

The appropriate ISCO solution is derived from 

= ro^ -6Mro + 8M^/\y^-3M\'^ ,{A10) 

where x = a/M. The solution of the above quartic func- 
tion is obtained as follows. For x > 0, we have 



r/sco = M {3 + V3x2 + A2 



(3 - A)(3 + A + 2 V3x' + A2) \ , (All) 



1/2 



and for x < 0, 

r/sco = Af {3 + ^3x2 + A2 



(3 - A)(3 + A + 2 v/3x2 + A2) \ , (A12) 



1/2 



where 



2U/3 



i + (i-x') 



(i + x)^/'V(i-x) 



1/3 



(A13) 



We can obtain the well known result, rjsco = A/ for 



a = M an d risco = 9Af for a = -Af from Eqs. (|Allj) 
and (A12|, respectively. 



2. Polar circular orbit 

For the polar orbit, we need to consider the limit 
2/ — >■ cx) in Eq. ( A8 1 . This means that we may solve 



the following equation. 



+3aVo^ -6a'^Mro + a^) . 



(A14) 



Here, we introduce two nondimensional variables, f 
and X as 



ro = Mxf, 
a = Mx- 



(A15) 
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Here we consider the case for x 7^ 0- Then, the equation 
to find the IS CO radius is written as 

= f^--6 — + 3f^+4 — 



the equation with the real coefficient, 



= (f^ - aif + 1) , 



(A19) 



+3f2-6- + l 



(A16) 



where we find 



The solutions are nice relations (we can find them from 
numerical method), and the solutions are given by 



1 



, exp(z6'i) , exp(— 16*1) , 



exp(i^2) , exp(~i6'2) 



(A17) 



where rg, 9i and 62 are real. The above solutions suggest 
that the 6th order equation can be reduced to 

= (f^ - aif + l)(f^ - + 1) 

x{f^ -a^f + l). (A18) 

Here, although ai is real, a2 and 013 are complex and 
complex conjugate each other. Therefore, we focus on 



ai = 2 



1/3 



+2 
2 



1/3 



X (1 - + «X\/2 - X^) 
(exp(i/3/3) + exp(-i^/3) + 1) ; 



tan/? = 



X 



Xv/2^ 

i-x^ 



(A20) 



As a result, we have the ISCO radius as the following. 



risco = m{ (1 - x' + *X V2 - x^) + (l - x' + ^X^a - x") + 1 



-1/3 



l-x2+*Xv/2-x2) +(1- 



1/3 



X^ + «xV2 - x^ 



-1/3 



We show the polar ISCO orbit for the case of a = 0.9M in Fig. [29l 



(A21) 
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